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Abstract 

Asymptotic behaviors are often of particular interest when analyzing 
Boolean networks that represent biological systems such as signal trans¬ 
duction or gene regulatory networks. Methods based on a generalization 
of the steady state notion, the so-called trap spaces, can be exploited to 
investigate attractor properties as well as for model reduction techniques. 
In this paper, we propose a novel optimization-based method for com¬ 
puting all minimal and maximal trap spaces and motivate their use. In 
particular, we add a new result yielding a lower bound for the number 
of cyclic attractors and illustrate the methods with a study of a MAPK 
pathway model. To test the efficiency and scalability of the method, we 
compare the performance of the ILP solver GUROBI with the ASP solver 
POTASSCO in a benchmark of random networks. 


1 Introduction 

Boolean network models have long since proved their worth in the context of 
modeling complex biological systems m- Of particular interest are number and 
size of attractors as well as their locations in state space since these properties 
often relate well to important biological behaviors. 

In this paper, we explore the notion of trap space, which is a subspace of 
state space that no path can leave, for model reduction and attractor analysis. 
After providing the relevant terminology, we demonstrate that trap spaces can 
be used for model reduction as well as to predict the number and locations of 
the system’s attractors. We then introduce the prime implicant graph as an 
object which captures the essential dynamical information required to compute 
trap spaces of a Boolean network. 

In practice, methods that exploit trap spaces can only be useful if their iden¬ 
tification scales efficiently with the size of the network. We provide an Integer 
Linear Programming (ILP) and an Answer Set Programming (ASP) formula¬ 
tion and present the results of a benchmark that compares the performances of 
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fl=Vi+ V2 
f2=Vi ■ V4 
f3=W -Vi 
f4 = V3 


/l(1101) = 1 + 1 = 1 
/2(1101) = 1-1 = 1 
/3(1101)=T-1 = 0 
/4(1101) = 0 = 1 


F(llOl) = 1101 
F(OOOO)= 0001 
F(OllO) = 1000 
F(llll) = 1100 


Figure 1: (left) An example Boolean network with four variables, (middle) 
An example for how the Boolean expressions are evaluated for a given state 
X = 1101. (right) Some examples for the image F{x) of a state x G S. 


the ILP solver GUROBI and the ASP solving collection potassco. Finally, we 
apply our methodology to a MAPK pathway model. 

This paper is an extended version of |12) . To allow for a more intuitive 
understanding of technical terms, we decided to reformulate some of the previ¬ 
ously presented notions in terms of the state space of Boolean networks. The 
extension consists of discussing both maximal and minimal trap spaces, an ASP 
formulation of the optimization problem, the benchmark and an extension of 
the MAPK case study. 

1.1 Background 

We consider variables from the Boolean domain B = {0,1} where 1 and 0 
represent the truth values true and false. A Boolean expression f over the 
variables V = {vi,..., n„} is defined by a formula over the grammar 

/ ::= 0 I 1 I n I 7 I /i • /2 I /i + /2 

where v G V signifies a variable, / the negation, fi ■ f 2 the conjunction and 
fi + /2 the (inclusive) disjunction of the expressions /,/i and / 2 . Given an 
assignment x : V —>■ B, an expression / can be evaluated to a value f{x) G B 
by substituting the values x{v) for the variables v G V. If f{x) = f{y) for all 
assignments x,y : V —>■ B, we say / is constant and write / = c, with c G B 
being the constant value. A Boolean network (V, F) consists of n variables 
V = {vi,... ,Vn} and n corresponding Boolean expressions F = {/i,.. •, /n} 
over V. In this context, an assignment a; : F —)■ B is also called a state of the 
network and the state space S consists of all possible 2” states. We specify 
states by a sequence of n values that correspond to the variables in the order 
given in V, i.e., x = 110 should be read as a:(ni) = l,x{v 2 ) = 1 and x{v 3 ) = 0. 
The expressions F can be thought of as a function F : S ^ S governing the 
network behavior. The image F{x) of a state x under F is defined to be the 
state y that satisfies y{vi) = fi{x). To illustrate these concepts we introduce a 
running example in Fig. 

The state transition graph of a Boolean network (V, F) is the directed graph 
(S', —>) where the transitions —> C S x S are obtained from F via a given 
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update rule. We mention two update rules here, the synchronous rule and its 
transition relation C S x S, and the asynchronous rule and its transition 
relation ^ C S x S. The former is defined hy x ^ y iS F{x) = y. To define ^ 
we need the Hamming distance A : S' x S' —>■ {1,..., n} between states which is 
given by A{x,y) := \{v €V \ x{v) ^ J/(u)}|. We define a: ^ y iff either x = y 
and F{x) = x or A{x,y) = 1 and A{y,F{x)) < A{x,F{x)). 

A path in (S, —>) is a sequence of states (xi ,..., a;i+i) with xt —>■ Xi+i, for 
1 < i <1. An non-empty set T C S is a trap set w.r.t. —> if for every x G T and 
y G S with X ^ y it holds that y G T. An inclusion-wise minimal trap set is 
also called an attractor of (S, —?►). Note that every trap set contains at least one 
minimal trap set and therefore at least one attractor. A variable u G H is steady 
in an attractor A C S iff x{v) = y{v) for all x,y G A and oscillating otherwise. 
We distinguish two types of attractors depending on their size. If A C ^ is an 
attractor and |A| = 1 then A is called a steady state and if |A| > 1 we call 
it a cyclic attractor. The cyclic attractors of {S, are, in general, different 
from the cyclic attractors of The steady states, however, are identical 

in both transition graphs because a; —> x iff F{x) = x in both cases. Hence, we 
may omit the update rule and denote the set of steady states by S'f. The state 
transition graphs of the running example is given in Fig. 


2 Methods 

2.1 Trap Spaces 

A subspace of S is characterized by its fixed and free variables. It may be 
specified by a mapping p : U where U CV is the subset of fixed variables, 
p{u) the value ot u G U and the remaining variables, V \ U, are said to be free. 
Subspaces are sometimes referred to as symbolic states or partial states. We 
specify subspaces like states but allow in addition the symbol to indicate that 
a variable is free, i.e., p = **10 means U = {x 3 ,X 4 } and piv^) = l,p{v 4 ) = 0. 
The set S* = S*{V) denotes all possible 3" subspaces. Note that states are a 
special case of subspace and that S C S* holds. We denote the fixed variables 
U of p G S* hy Up C V. A subspace p references the states 5'[p] := {x G S' | 
t/v G Up : x{v) = p{v)}. 

A trap space is a subspace that is also a trap set. Trap spaces are related to 
the seeds in [15] and stable motifs in m- Trap spaces are therefore trap sets 
with a particularly simple structure. They generalize the notion of steadiness 
from states to subspaces. We denote the trap spaces of (S, —>) by Sf^. Since 
every trap set contains at least one attractor, inclusion-wise minimal trap spaces 
can be used to predict the location of a particular attractor while maximal trap 
spaces may be useful in understanding the commitment of the system to a set 
of attractors. Hence, we define a partial order on S* based on whether the 
referenced subspaces are nested: p < q iS S'[p] C S[q]. The minimal and 
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maximal trap spaces are defined by 

min(5'^) :={p G S*_^ \ $ q & : q < p} 

max(5'^) :={p G \ $ q G : {q > p) A (^[g] ^ S')} 

where we use the condition S[g] 7 ^ S because otherwise S is a priori the (unique) 
maximal trap space of any network which is only informative in the special case 
that the network has a single attractor. The minimal and maximal trap spaces 
of the running example are illustrated in Fig. Note that in this case, for each 
attractor A C S of (S, ^) there is p G min(Si^) s.t. S[p] = A. 

2.2 Characterization of Trap Spaces 

Analogous to the characterization of steady states by the equation F{x) = x we 
show that trap spaces can be characterized by the inequality p > F[p\. 

To do so, we first dehne f[p] to be the expression obtained by substituting the 
values p{v) for the variables v G Up in /. The image F[p] of a subspace p under 
F = {/i,..., /„} is then the subspace q with Uq := {vi GV \ fi[p] is constant} 
and q{vi) := fi[p], for all Vi G Uq. The following theorem characterizes trap 
spaces. 

Theorem 1. A subspace p is a trap set of (S', —>■) if and only if p > F[p]. 

Proof, p G S* is a. trap set of {S, -a) iff there are no x G S[p] and y G S \ ^[p] 
s.t. X —>■ p. This is equivalent to /i(x) = p{vi) for all x G S'[p] and Vi G Up^ 
which is equivalent to p > T[p]. □ 

Note that the proof of Thm. is identical for the synchronous and asyn¬ 
chronous transition graphs. As a corollary we get that i.e, that trap 

spaces are (like steady states) independent of the update rule. We therefore 
write Sp instead of SA^. For our running example, note that the trap spaces of 
(5, '^) and {S, -^) are indeed identical, see Fig.[^ Note also that Sp A min(5'^) 
and so calculating all minimal trap spaces yields, in addition to the information 
on cyclic attractors, also all steady states. Example calculations are given in 
Fig. [31 

2.3 Applications 

2.3.1 Application 1: Model Reduction 

Let T C S' be a trap set. Denote by Sub(T) G S* the smallest subspace that 
contains T C S, i.e., Usub(T) '.= {v G V \ Vx,p G T : x{v) = y{v)} and 
Sub{T){v) := x{v) for x G T arbitrary. Note that, in general, the smallest 
subspace that contains T is a superset of T, i.e., S[Sm 6 (T')] A T. A natural 
model reduction technique is based on the observation that for any trap set 
T C S, the transitions of any path with an initial state Xi G T only depend on 
the reduced system {Vp, Fp) with p := Sub{T) and 

Vp — {v gV \ v ^ Up}, Fp := {/4p] \ fiG F, Vi G Vp}. 
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(b) The synchronous transition graph (S, —») 


Figure 2: The synchronous and asynchronous transition graphs of the run¬ 
ning example. The trap spaces (and the way they are nested) are given as a 
heat map. Minimal and maximal trap spaces are outlined by a black line and 
the attractors by a white line. Using the notation we get the following 
sets: S*. = 00^^, 1*0^, 1^01,1101}, min(,S^.) = {00^^, 1101} and 

max(S'l^) = {00**, 1***} and Sp = {1101}. 
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Figure 3: Three examples of how to compute the image of a subspace under F 
for the running example, p and q satisfy the the condition of Thm. and are 
therefore trap spaces. 


Intuitively speaking, the network (I^, Fp) is obtained by ’’dividing out” the fixed 
variables of p, see m for more details. Note that, by definition, 5 'm&(S'[p]) =p 
for any p G S*. In particular, trap spaces p have this property and so they 
naturally lend themselves for this reduction technique. 

2.3.2 Application 2: Cyclic Attractors 

The following theorem is based on the observation that a minimal trap space is 
either itself a steady state or it contains no steady state at all. 

Theorem 2. | min(S'J.)\5'i?| is a lower bound on the number of cyclic attractors 
of{S,^). 

Proof. Let p G min(5'^) \ Sp- Since S'[p] is a trap set, it contains an attractor 
A C 5'[p]. If A = {x} then x G Sp such that x < p, which contradicts the 
minimality of p. □ 

Furthermore, since Sub{S[p]) = p for p G Sp, we may conclude that some 
V G V \ Ugub{Slp]) must be involved in the cyclic behavior. In our running 
example, the subspace p = 00** is a minimal trap space and therefore contains 
only cyclic attractors in which either or va or both oscillate, see Fig. 

2.4 The Prime Implicant Graph 

In this section we propose a method for computing trap spaces. The idea is to 
translate the task into a hypergraph problem in which trap spaces are repre¬ 
sented by a sets of arcs that satisfy certain constraints. We consider a directed 
hypergraph in which each arc corresponds to a minimal size implicant of fi or 
fi, for some Vi G V. Minimal size implicants are also called prime implicants, 
see e.g. [T]. Here, we define the following slight variation: For c S B, a c- 
prime implicant of a non-constant expression / is a subspace p G S* satisfying 
f[p] = c, and f[q] ^ c for all q > p. For a constant fi = ewe define that p with 
Up := {ui} and p{vi) = c is its single prime implicant. In the running example. 
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11** satisfies /i[ll**] = 1. But it is not a 1-prime implicant of /i, because 
1 *** > 11** and [1***] — 1- 

The set of all prime implicants of a Boolean network is denoted hy V = 
V{V,F) and consists of all {p,c,Vi) S S'* x B x 1/ such that p is a c-prime 
implicant of fi. The prime implicant graph is the directed hypergraph (A/", „4) 
where Af = Af{V) := {p G S* | \Up\ = 1} consists of all subspaces with \Up\ = 
1 (corresponding to literals in propositional logic). To define the arcs A = 
A{V,F) C2-^ X 2-^ we observe that a subspace p can be written (uniquely) as 
the intersection between k := \Up\ subspaces pi,... ,pk such that \Up.\ = 1 for 
all i. 

The arcs A = A{V, F) C 2^ x 2^ are defined by the mapping 
a:V ^2^ x2^, {p,c,Vi) ^ ({pi, • ■ ■ ,Pfc}, {?}) 


where 

(1) pi,... ,pfc is the decomposition of p into k = \Up \ subspaces pi with \ Up. | = 

1 . 

( 2 ) q is defined by Uq := {vi} and q{vi) := c. 

The prime implicant graph has exactly one arc for every prime implicant, i.e., 
A := {q;(p, c, Vi) I (p, c, Vi) G V}. The head of an arc a = ({pi,. .. ,Pfe}, {(?}) is 
denoted by H{a) := q, and its tail by T{a) := p where p is the intersection of 
Pi,... ,Pfe. The prime implicant graph of the running example is given in Fig.[^ 

2.5 Prime Implicants and Trap Spaces 

Now we establish a relationship between subsets of A and trap spaces. To 
do so we need the notions of consistency and stability. A subset A C A is 
consistent if for all oi, 02 G A and v G UH(ai) H UH(a 2 ) it holds that F[{ai)(v) = 
H{a2){v). If A = {ai, ..., 0 ^} C A is consistent then the intersection between 
7J(ai),..., F[{am) is non-empty, called the induced subspace of A and denoted 
by H{A). For the special case A = 0 we define H{A) by Uh{a) ■= A. 
subset A C A is stable if for every a £ A there is a consistent subset Ba 'F A 
such that T{a) > F[{Ba). Intuitively, in this case the requirement T{a) for each 
implication a G A to become effective is met by some assumptions Ba- The 
stable and consistent subsets of A for the running example and their induced 
subspaces are given in Fig. The central idea for the computation of S*p is 
given in the next result: 

Theorem 3. p G S* is a trap space if and only if there is a stable and consistent 
A C A such that H{A) = p. 

Proof. The statement of the theorem is true ii Up = 0 because A = 0 is stable 
and consistent and F[ (A) = p by definition. 

Hence, assume Up ^ 0. Let Vi G Up. Since p > F[p] it follows that fi[p] = 
p{vi) and hence that there is a. p{vi)-pTime implicant qi of ft that satisfies qi > p. 
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1, I—>■ di 
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^ = {fll} 

B = {01,08} 
c = (oi, 08 , Oioj 
H(Di) = 1101 


ff(A) = 1*** 
H{B) = 1*0* 
H{C) = 1*01 
H{D2 = 1101 

(c) 


Di = {01,04,08,010} 

D2 = {02,04,08,010} 

D3 = {01,02,04,08,010} 
HiDs) = 1101 


Figure 4: (a) The complete set V of prime implicants of the example network, 
given in (p, c, Vi) notation, (b) The prime implicant graph (Af, .4) of the running 
example. The eight nodes are grouped into four pairs that belong to the same 
variable Vi. Gray discs represent the groups and are not elements of Af. Nodes 
that correspond to positive literals are drawn in black and negated literals in 
white. Hyperarcs are represented by several arcs having a common arrowhead. 
The colors indicate three of the five stable and consistent arc sets, (c) The 
non-empty stable and consistent arc sets and their induced trap spaces. 
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The set A := {a(qi,p(vi),Vi) | Vi G Up} is, by definition, consistent and satisfies 
H{A) = p. But it is also stable because T(a) > p for all a G A. Let 0 ^ A C A 
be stable and consistent. Then there exists a G A s.t. H{a) = H{A){vi) for all 
Vi G Uh{a)- Hence H{A) > F[H{A)] and H{A) G S*p. □ 

The following corollary points out that extremal (i.e., minimal or maximal) 
stable and consistent arc sets correspond to extremal trap spaces. As a partial 
order on the subsets of A we use set-inclusion. 

Corollary 1. (1) Ifp G max(Sp) then there is a non-empty minimal stable and 
consistent A C A s.t. H{A) = p. (2) If p G min(S'^) then there is a maximal 
stable and consistent A G- A s.t. H{A) = p. 

Note that the inverse relationship (maximal arc sets induce minimal trap 
spaces and the other way around) stems from the fact that the order on S* con¬ 
siders the free variables while the order on arc sets is in terms of fixed variables. 
Note that stable and consistent arc sets are a generalization of the so-called 
self-freezing circuits which were described and studied in [3l [10]. These cir¬ 
cuits are based on canalizing effects of F which correspond to prime implicants 
{p,c,Vi) G V that satisfy \Up\ = 1. Self-freezing circuits therefore correspond to 
trap spaces whose stable and consistent AC A satisfy \UT{a) \ = 1 for all a G A. 

2.6 Computation of Trap Spaces 

In this section we formulate a 0-1 optimization problem to compute all extremal 
stable and consistent A C A and therefore all extremal trap spaces. To solve it 
in practice, we suggest using solvers for Integer Linear Programming (ILP) or 
Answer Set Programming (ASP) and give a reformulation of the constraints in 
terms of each language. 

As a preliminary step the set 7^(H, F) has to be enumerated, see e. g. [9|. Al¬ 
though the number of prime implicants may grow exponentially with the number 
of variables that the expression depends on, see e. g. [T], we found that for typ¬ 
ical biological models these dependencies are so small that the enumeration of 
V is negligible compared with finding consistent and stable arc sets. 

We now formulate the 0-1 optimization problem. For every arc a = 
{p,c,Vi) G V we introduce a variable Xa G {0,1} that indicates whether or 
not o is a member of the set AC A that we want to compute. We denote these 
variables by A •.= {xa \ a GV}. In addition, we introduce for every Vi G V two 
variables i/°, yj G (0,1} that indicate whether Vi is fixed in the trap space and, 
if so, what value it takes. We denote them by Y := {yf | c G 18 , 1 ;^ G V}. For 
any c G B, we require = 1 if and only if Vi G Uh{a) and H{A){vi) = c. To 
encode this requirement, we use the logical constraints 

Vi V c G B,r)j G y, (Cl) 

aeBf 

where Bf := {a G A \ {I'i} = UH{a),II{a){vi) = c} denotes the arcs inducing 
Vi to take the value c, and ^ and V are the standard logical connectives for 
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implication and disjunction. To enforce that the set A := {a G | Xa = 1} is 
stable and consistent we add the following constraints (C2) resp. (C3): 


Xa ^ for all a G G 17r(a), (C2) 

Vi V y], for all Vi G V. (C3) 


The maximal resp. minimal stable and consistent A A correspond to all 
solutions of 

max ^ Xa, s.t. (C1),(C2),(C3), (0-1 max) 

respectively 

min ^ Xa, s.t. (C1),(C2),(C3) and 1 < ^ Xa, (0-1 min) 

Xa&X Xa&X 


where the additional constraint for the minimal solutions forbids the empty set. 


2.6.1 ILP Formulation 

We now reformulate the above constraints by linearizing the logical operators. 
The resulting problem can be solved with standard ILP solvers. The first con¬ 
straints, (ILPl), enforce that j/f if there is no arc targeting it: 

yi < '^ Xa and Va G : Xa < y^, (ILPl) 

for every c G B and Vi G V. The next constraints enforce stability and consis¬ 
tency, respectively: 

for all a G .4, Ui G t/T(a), (ILP2) 

y^ + < 1, for all Vi G V. (ILP3) 

A PYTHON implementation using GUROBI [5] is available from m- Since ILP 
solvers usually do not return multiple solutions we suggest to iteratively add no- 
good-cuts that make the current solution infeasible but do not otherwise affect 
the feasible region, until there are no more solutions. Suppose s : XLlY —>■ {0,1} 
is a solution. Cuts for the maximization and minimization of (0-1), respectively, 
are 


1 < Xa, 

XaGG 

^ Xa< \E\ - 1, 

Xa€:E 


where G := {xa G X \ s{xa) = 0} 
where E := {xa G X \ s{xa) = !}■ 
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2.6.2 ASP Formulation 


We now reformulate the constraints (C1)-(C3) as an answer set program, see e.g. 
[^. To encode the arcs A we introduce two ternary predicates, head(v,c,ID) 
and the tail(v,c,ID), where v refers to some w € F, c to a value in c G B 
and ID is an index that determines to which arc a tail or a head belongs. Each 
a G A is then translated into a number of so-called facts by stating all the 
tail elements and the unique head element it consists of. For example, an arc 
03 = (00**, 0, ui) of the running example in Fig. [^becomes 

head(vl,0,a3). tail(vl,0,a3). tail(v2,0,a3). 

Note that the index a3 links the data together. In the generate-and-test fashion 
of defining ASP problems we generate all possible subsets of A and introduce 
an unary predicate x(ID) to indicate whether the arc with index ID belongs to 
the solution A C A or not. It encodes the variables G A in the formulation 
of the (0-1) problem above. 

{x(ID) : head(v,c,ID)}. 

The ASP formulation does not require the auxiliary variables Y and hence can 
do without (Cl). The stability (C2) is translated into the hlter (ASP2) which 
forbids the existence of an arc (with identifier IDl) such that one of its tail nodes 
is not also the head of another arc (with identifier ID2). The consistency (C3) is 
translated into the filter (ASP3). It forbids that two arcs (with identifiers IDl 
and ID2) target the same v but at different values 0 and 1. 


x(IDl), tail(v,c,IDl), not x(ID2): head(v,c,ID2). (ASP2) 


x(IDl), x(ID2), head(v,l,IDl), head(v,0,ID2). (ASP3) 

To compute multiple solutions is built into ASP solvers and the solving 
collection potassco also features the option to find set-inclusion minimal or 
maximal solutions with respect to the predicates that are true. For the problem 
at hand we optimize over the predicate x(ID). To forbid the empty solution 
when minimizing we use the additional filter ” ; - {x(_)} 0. ”. A python 
implementation using the solving collection POTASSCO is available at m- 

2.6.3 Benchmark 

To test the efficiency and scalability of the ASP and ILP formulations we cre¬ 
ated random Boolean networks and recorded the time necessary to compute 
min(S'J.) and max(S'^) with each solver. For the results to be reproducible 
we decided to use the function generateRandoniNKNetwork of the R package 
BOOLNET [13]. It takes the parameters n and k which specify the number of 
variables n = \V\ and number of variables k that each f G F depends on. 
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POTASSCO ASP 


GUROBI ILP 



(a) Running times for the different solvers 




(b) Number and sizes of trap spaces 

Figure 5: (a) The average running times of POTASSCO and GUROBI for comput¬ 

ing min(S'^) and max(S'^). Each plot also contains the total number of prime 
implicants \V\. Discontinued plots indicate a time-out of the respective solver, 
(b) The plot on the left shows how many variables \Up\ are (on average) fixed in 
a minimal and maximal trap space p. The dashed line is a plot of f{n) = n and 
therefore indicates how close a trap spaces is to being a steady state. The plot 
on the right shows the average number of minimal and maximal trap spaces. 
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In addition to n and k there are different configurations for the topology, the 
linkage and the bias in the generation of the truth tables of the Boolean ex¬ 
pressions. We decided to generate networks whose in-degree follows a Poisson 
distribution (topology="homogeneous") with mean and variance equal to fc = 3 
and left the other parameters at their default values (linkage="unif orm" and 
functionGeneration="uniform"). This setup allows variables with high in¬ 
degrees (hubs) as well as low in-degrees (e.g. inputs, cascades, outputs) that 
frequently appear in models of biological networks. 

For each n, starting from n = 50, we generated 50 networks and called 
GUROBI and POTASSCO to find all minimal and maximal trap spaces of each. For 
each network we recorded the number of prime implicants, number of minimal 
and maximal trap spaces, average number of fixed variables in the trap spaces 
and the average time for each solver to compute them. 

We treated each solver and whether minimal or maximal trap spaces are 
computed, as an independent computation and allowed a time limit of 10 min¬ 
utes for each. When a computation failed (time-out or out of memory) 25 times 
for the same n we removed it from the benchmark loop. To solve the ASP 
problems, we used GRINGO version 3.0.5 and CLASP version 3.1.1 with the con¬ 
figuration —dom-pref=32, —heu=domain, —doni-mod=6 (subset minimality) or 
—dom-mod=7 (subset maximality). To solve the ILP problems we used used the 
PYTHON interface to GUROBI, version 5.0. The executions were ran on a Linux 
desktop PC with 30 GB RAM and eight CPUs with 3.00GHz. The results are 
given in Fig. 

The first observation is that the POTASSCO formulation performs better than 
GUROBI formulation for both maximal and minimal trap spaces. The time 
required to compute maximal trap spaces appears to grow linearly while the 
one for minimal trap spaces grows exponentially. A reason for why our ILP 
formulation is less efficient than the ASP formulation might be that the ILP 
solutions are constructed iteratively, by forbidding the last found solution, while 
the ASP solutions are computed in a single execution of CLASP. 

The two plots at the bottom of Fig. give some statistical information 
about the number of trap spaces and their sizes in terms of fixed variables. The 
number of fixed variables in a maximal or minimal trap space appears to be a 
fixed fraction of n and the number of maximal trap spaces is constant while the 
number of minimal trap spaces grows linearly with n. 


3 Application to a MAPK pathway model 

We computed the extremal trap spaces for a network that models the influence of 
the MAPK pathway on cancer cell fate decisions, as published in [7]. It consists 
of 53 variables that represent signaling proteins, genes and phenomenological 
components like proliferation or apoptosis. We found that there are 18 minimal 
trap spaces, 12 of which are steady states. Hence, following Application 2 
in Sec. |2.3[ there are at least six cyclic attractors whose properties can be 
comprehensively investigated using the six corresponding reduced models. The 
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Figure 6: (a) The six trap spaces min(5'p) \Sf = {pi, ■ ■ ■ jPe} in terms of fixed 
variables in black and white and free variables in red. (b) We mapped the fixed 
and free variables of p 4 onto the interaction graph of the MAPK network. The 
corresponding reduced network consists of \V\ — \UpJ = 53 — 45 = 8 variables 
and is given on the right. 


trap spaces min(S'^) \ Sp as well as one reduced model are given in Fig. Next 
we ask how many attractors A C S'[p] there are for each p S min(S'^) \ Sp- We 
used BNS [5] to compute all attractors of (5,-») and GENYSIS [3] to compute 
the attractors of {S, ^). 

It turns out that for the asynchronous update, each trap space p S min(5'^) 
contains exactly one attractor. In addition, the attractor A C S'[p] satisfies 
Sub{A) = p for each p, i.e., the fixed and free variables of p correspond exactly 
to the variables that are steady and oscillating in its attractor A. We obtained 
this result by computing the attractors of the reduced networks (V^, Fp) of each 
p G min(5'^) with GENYSIS. The question whether there is an attractor outside 
of Upgmin(s* )5'[p] Can not be answered with this method because the unreduced 
MAPK network is too large for GENYSIS. 

All attractors of the synchronous transition graph can be computed with 
BNS, even for the unreduced MAPK network. We found that there are 28 cyclic 
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Table 1: For each of the 9 maximal trap spaces we counted how many steady 
states and cyclic attractors in the synchronous and asynchronous transition 
graphs they contain. The first eight columns correspond to the inputs at either 
0 or 1 and the last column refers to the maximal trap space in which PI3K and 
GABl are both steady at 1. Note that {TGFBR_stimulus = 1) alone ensures 
that the system will end up in a steady state (zero cyclic attractors reachable) 
and that most steady states (10 out of 12) are inside the PI3K, GABl trap 
space. 


maximal trap spaces 

EGF-stim 

FGF-stim 

TGF-stim 

DNA-dmg 

PI3k, G 

0 

1 

0 

1 

0 

1 

0 

1 

1 

Steady states 

8 

4 

8 

4 

4 

8 

6 

6 

10 

Cyclic attractors in (S',-») 

24 

4 

17 

11 

28 

0 

22 

6 

19 

Cyclic attractors in (S,^) 

2 

4 

2 

4 

6 

0 

3 

3 

6 


attractors. Contrary to the asynchronous case, some minimal trap spaces do 
contain more than one attractor and Sub{A) < p does frequently hold. For 
example, pi contains four cyclic attractors and although |17pJ = 12 between 19 
and 24 variables are steady in the attractors. In addition, 18 attractor are not 
contained in any minimal trap space. 

We then computed the maximal trap spaces and found that | max(5'p)| = 9. 
The MAPK network has four input variables v G V, namely EGFRstimulus, 
FGFR3stimulus, TGFBRstimulus and DNA^damage, that satisfy /„ = v. 
Each input generates two maximal trap spaces p^,p^ defined by Upc = {u} and 
p'^(u) = c, for c S B. Eight of the maximal trap spaces are therefore explained 
by the inputs. The additional q G max(S'^) is defined by Uq := {PI3K, GABl} 
and q{PI3K) = q{GABl ) = 1. A summary of the commitments of the different 
maximal trap spaces to steady states and cyclic attractors is given in Tab. 

4 Discussion 

In this paper we propose to use trap spaces for the prediction of a network’s 
attractors and for model reduction. We propose a novel, optimization-based 
method for computing trap spaces that uses ILP or ASP solvers. Its input is 
the prime implicant graph rather than the full state space (which is necessarily 
exponential in E). The method can be extended from Boolean to multi-valued 
networks by generalizing the notion of prime implicants from Boolean to multi¬ 
valued expressions. Thm. the characterization of trap spaces, holds not only 
for synchronous or asynchronous but for any update rule (see [6j for other 
update rules) and also stochastic simulations. Note also that the prime implicant 
graph is similar to the process hitting graphs in |14j . certain prime implicants 
correspond for example to actions, but different in that it attempts to capture 
a single model exactly rather than to approximate a family of models. We are 
not aware that process hitting has been used to compute trap spaces. 
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The benchmark demonstrates that the prime implicant method for com¬ 
puting trap spaces scales well with the number of variables of a network. All 
maximal and minimal trap spaces, including steady states, of a network with 
between 300 and 400 variables and poisson-distributed connectivity with fc = 3, 
are computed with an average running time on the order of minutes. A com¬ 
parison with a brute force approach for computing trap spaces, by enumerating 
subspaces, and the circuit-enumeration approach presented in m would be 
desirable to further assess the efficiency of this approach. 

The results for the MAPK network suggest that the usefulness of minimal 
trap spaces in approximating attractors depends a lot on the update rule. For 
the asynchronous update, we observe that (1) eachp G min(5'^) contains exactly 
one attractor A and that (2) Up corresponds exactly to the steady variables of A. 
For the synchronous update, none of these properties are met and, additionally, 
there are 18 attractors that lie outside of any minimal trap space. It appears 
to us that the trap space method for predicting attractors is, in general, less 
precise for synchronous transition graphs. 

Regarding the use of maximal trap spaces it is worth noting that although 
they usually contain several attractors, they are likely to play an important role 
in decision making and processes like cell differentiation. An example is the 
observation that TGFBR = 1 eliminates the possibility of sustained oscillations 
and guarantees that the system ends up in a steady state. 

Currently, we are working on a method that combines computing trap spaces 
with reduction techniques and model checking to decide, for a given network, 
whether properties ( 1 ) and ( 2 ) hold. 
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